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Abstract 



We present a high order perturbative computation of the renormahzation con- 
stants Zy, Za and of the ratio Zp/Zs for Wilson fermions. The computational 
setup is the one provided by the RI'-MOM scheme. Three- and four-loop expan- 
sions are made possible by Numerical Stochastic Perturbation Theory. Results 
are given for various numbers of flavors and/or (within a finite accuracy) for 
generic rij up to three loops. For the case = 2 we also present four-loop re- 
sults. Finite size effects are well under control and the continuum limit is taken 
by means of hypercuhic symmetric Taylor expansions. The main indetermination 
comes from truncation errors, which should be assessed in connection with con- 
vergence properties of the series. The latter is best discussed in the framework 
of Boosted Perturbation Theory, whose impact we try to assess carefully. Final 
results and their uncertainties show that high-loop perturbative computations of 
Lattice QCD Renormalization Constants (RC's) are feasible and should not be 
viewed as a second choice. As a by-product, we discuss the perturbative expan- 
sion for the critical mass, also for which results are for generic ri/ up to three 
loops, while a four-loop result is obtained for ri/ = 2. 



1 Introduction 



Lattice Perturbation Theory (LPT) has been for a long time the only available tool for the 
computation of Lattice QCD Renormalization Constants (RC's). By now, non-perturbative 
computations are preferred. We should stress, however, that there is no theoretical obstacle 
to the perturbative computation of either finite or logarithmically divergent RC's like (for 
example, those for quark bilincars or their ratios). The main difficulties are of practical 
nature. The first one is that LPT is technically very hard, much harder than Perturbation 
Theory (PT) on the continuum [1]. Therefore, computations arc often performed only at one 
loop. This is a serious limitation, which is made even more severe by the bad convergence 
properties of LPT. To take care of this problem Boosted Perturbation Theory (BPT) and/or 
of Tadpole-Improved Perturbation Theory (TIPT) [2] is often used. There is quite a consen- 
sus on the fact that, at one loop the impact of BPT (and/or TIPT) is often important. On 
the other side, there is no clear-cut result on the actual control on these procedures. One 
should always keep in mind that convergence properties of the series are the real issue and 
assessing them from a one-loop computation is of course impossible. Other improvement 
schemes have been in recent years proposed, which aim at rcsumming some leading contri- 
butions [3]. A different approach to the computation of RC's in Lattice QCD is a completely 
non-perturbative one. In this case, one needs an intermediate scheme, which is eventually 
matched to the MS scheme (the one in which phenomenologists are most interested in) by 
a continuum perturbative computation. Popular intermediate schemes are the Regulariza- 
tion Independent (RI'-MOM) [4] scheme and the Schrodinger Functional (SF) [5] scheme. A 
non-perturbative computation eliminates the truncation errors. On the other side, one needs 
to face all the difficulties inherent in numerics. Among these, the high computational effort 
of unquenched (or partially quenched) Lattice QCD simulations. In practice it is sometimes 
extremely hard to get a good signal for realistic simulation parameters. For this reason LPT 
is still necessary, either for comparison or because the only feasible approach. 

In recent years the technique of Numerical Stochastic Perturbation Theory (NSPT) has 
been introduced (for an extended introduction - which in particular covers the unquenched 
version - see [6]). NSPT is a numerical implementation of Stochastic PT [7]. It is a nu- 
merical tool which enables to perform LPT computations with no reference whatsoever to 
diagrammatics. By making use of NSPT we can compute Lattice QCD RC's to high orders, 
which here means 3 (or even 4) loops. At these orders the use of BPT enables to assess the 
convergence properties of the series and better control the truncation errors. Since we neces- 
sarily work on a finite lattice, finite-volume and scaling violation effects have to be assessed 
carefully: this can be done. A careful extraction of continuum limit is one of the good point 
of the approach: the solution comes from what we call hypercubic symmetric Taylor expan- 
sions. Another nice feature comes from the fact that we can work directly in the massless 
hmit (which is also where RC's are usually defined). This also eliminates the need of an 
expensive chiral extrapolations. Finally, perturbative computations offer the possibility of a 
stronger analytical control, knowing the dependence on the coupling and on the number of 
fiavors. 

The main message of this paper is that high-loop perturbative computations of Lattice 
QCD RC's are feasible and should not be seen as a second choice. In particular, having 
both the perturbative and the non-perturbative determinations of a RC's gives a valuable 
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comparison. This is not at all academic. As a matter of fact, non-perturbative determinations 
are based on assumptions that are only proved in PT. 

This is the first of a couple of papers which deal with the NSPT perturbative compu- 
tation of Wilson quark bilinears. Renormalization conditions are fixed by the RI'-MOM 
prescriptions. Here we can make comparison with a non-perturbative determination [8].^ In 
this paper we will concentrate on the determination of finite RC's: Zy, Za and the ratios 
Zp/Zs and Zv/Za ^ for (unimproved) Wilson fermions. As a by-product, we also obtain the 
expansion for the critical mass. Results are given for various numbers of fiavors. At three 
loops some results are even given (to a finite accuracy) for generic nj. Instead, we present 
fourth loop results for the Uf = 2 case only. A different paper will deal with the compu- 
tation of logarithmically divergent RC's for quark bilinears (in particular, the RC for the 
scalar current Zs — Z~^, which is phenomenologically relevant for the determination of the 
quark masses). This deserves some extra caution, since dealing with anomalous dimensions 
requires a peculiar care for finite volume effects. 

The paper is organized as follows. In section 2 we recall the basic definitions of the 
renormalization scheme to which we adhere, while in section 3 we discuss some technical 
details of our computations. Section 4 introduces the main tool which is needed to extract 
the continuum limit (the already mentioned hypercubic symmetric Taylor expansions): this is 
done by discussing the (prototype) computation of the quark propagator. Section 5 contains 
our results: first we discuss the finite ratios Zp/Zs and Zv/Za, for which we can fit three 
loops results for generic Uf] then we move to Zy and Za (results are given at three loops for 
Uf = and at four loops for Uf = 2); finally we present a by-product of our computations, 
i.e. the critical mass to three loops (again, actually four in the case Uf = 2). In section 6 
we discuss the general features of computations dealing with an anomalous dimension (this 
sets the stage for what will be discussed in a following paper [9]). In section 7 we deal with 
resummations and convergence properties of our series and finally section 8 contains our 
conclusions and perspectives for future applications. 



2 The RI'-MOM renormalization scheme 

In order to compute renormalization constants we adhere to the RI'-MOM scheme. This 
is one of the so-called physical schemes^ (as opposed to the more popular MS scheme) and 
goes back to the MOM scheme of [10]. It became very popular after the introduction of 
non-perturbative renormalization in [4] . RI emphasizes the regulator independent nature of 
the scheme, which in particular makes the lattice a viable regulator. The prime denotes a 
renormalization condition for the quark field which is slightly different from the original one. 
All the details on this scheme can be found, for example, in [11]. In the following we only 
introduce the definitions which are relevant for our application. 

The basic quantities of our computation are of quark bilinears between external quark 

^The comparison will be made for given values of the coupling (/? = 5.8) and number of flavors (n/ = 2). 

^We will explain below why the computation of Zy, Za and Zv/Za are not tautological here. 

^One should nevertheless keep in mind that the name physical is actually misleading in the case of QCD. 
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states at fixed (off-shell) momentum p 

dx {p\ ^jj{x)Tilj{x) \p) = Gr{pa). (1) 



Here F stands for any of the 16 matrices, which provide the standard basis of the Dirac space 
(Dirac ind ices will often be suppressed). We adopt the usual naming convention for the 
bilinear: the scalar {S) is defined by F = 1, the vector {V) by F = 7^, the pseudoscalar (P) 
by F = 75, the axial (A) by F = 7^75 and the tensor (T) by F = a^^, = 1/2 [7^, 7i/]- Above, 
we made explicit the dependence on the lattice spacing a, which serves as a regulator. Later 
we will use the notation p — pa. 

Being these quantities gauge-dependent, a choice for the gauge condition has to be made. 
We will focus on computations in the Landau gauge. Prom a numerical point of view, this 
gauge condition is easy to fix on the lattice. On top of that, one does not need to discuss the 
gauge parameter renormalization. It also gives some extra bonus: the anomalous dimension 
for the quark field is zero at one loop. 

We can trade the Gr{pa) for the amputated function Fr(pa) {S{pa) is the quark propagator) 

Gripa) Tripa) = S-\pa) Gr{pa) S-\pa). (2) 

The Tripa) are eventually projected on the tree- level structure by a suitable operator Por 

Oripa) = IV (Po, Tripa)) . (3) 

Renormalization conditions are now given in terms of the Or (pa) according to 

Zoril^a,gia)) Z'^i/ia, gia)) Oripa) 2_ ^ = 1- (4) 

p — M 

Here the Z's depend on the scale via the dimensionless quantity fia, while the dependence 
on gia) will be expanded in FT. One should keep in mind from the very beginning that we 
will be eventually interested in the a — > limit of the Z's. The quark field renormalization 
constant Zg, which enters the above formula, is defined by 



7 r r w -1 Tr(^5 \pa)) 
Zgiiia^gia)) = -i- (5) 



The original RI scheme (without a prime), would have a derivative with respect to p^, 
instead. 

In order to get a mass-independent renormalization scheme, one imposes renormalization 
conditions on massless quarks. In Perturbation Theory this implies the knowledge of the 
relevant counterterms, i.e. the values of the various orders of the Wilson fermions critical 
mass. One- and two- loop results are known from the literature [12]. Third (and fourth) loop 
have been computed by us as a (necessary) by-product of the current computations: results 
are reported in section 5 (the three loop result in the Uf = 2 case has already been reported in 
[6]). Notice that the situation in the non-pcrturbative framework is more cumbersome with 
respect to staying in the massless limit. The determination of the critical mass is in a sense 
the prototype non-perturbative computation of a(n additive) renormalization constant. As 
it is well known, this is also a matter of principle: being a power-divergent renormalization. 
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the critical mass itself can not be computed in Perturbation Theory in the continuum limit. 
Still, from a numerical point of view the massless limit is always reached by an extrapolation 
procedure, which is usually a major source of error in non-perturbative determinations of 
RC's for Lattice QCD. 

A great advantage of working in the RI'-MOM scheme is that the relevant anomalous 
dimensions are known to 3 loops [11]. One is usually ready to admit that getting the 
logarithms is the easy part in the computation of a renormalization constant, while fixing 
the finite parts is the hard part of the work. As we will see, the situation is, to a certain 
extent, the opposite in the case of NSPT. We actually take for granted the logarithms 
(they are fixed by the choice of the scheme) and mainly concentrate on the computation of 
finite parts. As it is discussed in section 6, finite size effects open anyway the backdoor to 
corrections to the logarithmic contributions. Being 3 loops the order to which anomalous 
dimensions are known, this is also the order at which we can push our computations for 
every observable which has a non-vanishing anomalous dimension. On the other side, the 
finite RC's we will be concerned with in the present paper are in principle not constrained 
by anything but numerical precision, and that is why we pushed the computation of these 
quantities to an even higher order (4 loops, at the moment). 



3 Some technical details of our computations 

The lattice formulation we use in this work is defined by the plain Wilson action for gauge 
fields and plain {i.e. unimproved) Wilson fermions. As said, our computational tool is NSPT 
[6]. Here we only point out those technical details which are relevant to the present computa- 
tion. In its actual implementation, NSPT shares a few ingredients which are common to any 
lattice simulation. The main peculiarity is the representation of the fields as an expansion 
in the coupling constant, i.e. 

n 

u,{x)^i+j:(3-'^uj:\x). (6) 

i=l 

As it is apparent from the formula above, our preferred expansion parameter is the inverse 
of the lattice parameter (3 = 2Nc/go, Nq being the number of colors and go = g{a) the 
bare lattice coupling; thus /3^2 is proportional to qq. In our case a three-loop computation 
requires n — 6, while for four loops (which is at the moment the maximum order for which 
we report results in the Uf — 2 case) one needs n — 8. The proliferation of fields results 
in the request of a bigger amount of memory than in ordinary (non-perturbative) lattice 
QCD simulations. It is of course relevant also in terms of computing power: the algorithm 
is dominated by order- by-order multiplications, i.e. the number of floating-point operations 
grows as n(n — 1)/2. While this could appear as a big overhead with respect to ordinary non- 
perturbative dynamics, this is actually not true. In particular, in unquenched NSPT (like 
in any fermionic simulation) the basic building block is the inversion of the Dirac matrix, 
for which the perturbative nature of the computation results in a closed recursive algorithm, 
which is fairly well implemented (see [6]). On top of that, as we will discuss a bit more 
later, there is no need for an extrapolation to the chiral limit. As a result, NSPT fermionic 
computations are actually less demanding than non-perturbative counterparts. 
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As in many non-perturbative numerical computations, it is worth producing fairly decor- 
related configurations and store them for different subsequent measurements. A 32^ lattice 
(both at three and at four loops) fits well on an APEmille crate. At the same orders, a 16^ 
lattice can be managed by small PC-clusters or even by a robust (but nowadays standard) 
PC. While the first case is treated by our TAO"^ codes, the second is implemented in the 
framework of a by now well-established C"*""*" NSPT package. 

The number of flavors Uf enters the computations as a parameter, i.e. one has to perform 
different simulations for different n/'s. In Perturbation Theory each order has a trivial 
polynomial dependence on n/, so that one can fit the nf dependence. Uf = has by now 
been simulated both on 16*^ and on 32^ lattices: results have been used to assess finite-size 
effects. The unquenched cases have been simulated on the bigger (32"^) lattice: Uf — 2 is 
the case for which we have the largest number of configurations, while we also have several 
tenths of configurations for both Uf — 3 and n/ = 4. As a result, at the moment we are not 
going to quote every result for any Uf. In particular, four-loop results are at the moment 
only given in the case n/ = 2, for a reason that will be clear in a moment. 

As already stated, one good feature of NSPT computations is the fact that one can stay 
at the chiral limit. As we have already pointed out, the computation of the Wilson fermion 
critical mass was in a sense the prototype computation of a non-perturbative (additive, in 
this case) renormalization constant. It is also the prototype of a power-divergent renormal- 
ization, which can not be safely computed in PT in the continuum limit. On the other 
hand, no numerical simulation can be performed at kcriucai (we adhere to the common non- 
perturbative notation of quoting the hopping parameter rather than the mass of the quark): 
the chiral limit is always reached by means of a convenient chiral extrapolation. 

In Perturbation Theory one corrects for the additive quark mass renormalization by 
plugging in critical mass counterterms order by order. This is exactly what we do in NSPT. 
We were ready to start our simulations straightaway at three-loop order, which requires the 
knowledge of the critical mass up to two loops, and this is exactly what can be taken from 
the literature [12]. Each subsequent order asks for an iterative procedure: one computes 
the critical mass at the n*^ order (from n*^-ordcr simulations) and then plug it in the (n + 
l)*^-order simulations. In particular, for the case n/ = 2 our determination of the three- 
loop critical mass was good enough to plug it into four-loop simulations. The statistics we 
collected for the other values of are at the moment not sufficient to safely aim at the same 
accuracy. 

The RI'-MOM scheme renormahzation has been discussed in a generic covariant gauge 
[11]. We have already stated that our computations were performed in Landau gauge and 
stressed what the advantages of a such a choice are. From the point of view of computer 
simulations fixing the gauge to Landau in NSPT simply requires the ordcr-by-order imple- 
mentation of a well-known (FFT-accelerated) iterative procedure (for details, see [6]). It is 
worth stressing that in the NSPT framework also a peculiar implementation of the Faddeev- 
Popov mechanism is possible (see [13] for an application): by the same trick which enables 
us to treat the fermionic determinant we can manage the Faddeev- Popov determinant, with- 
out the inclusion of ghost fields. Still, we can compute in any covariant gauge with gauge 

'^TAO is the APE-dedica.ted programming language. 
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parameter ^ 7^ 0, i.e. Landau gauge is the only one which is not viable (apart from an 
extrapolation procedure). While the generic covariant gauge NSPT simulation has a (mod- 
erate) computational overhead, Landau gauge-fixing has a delicate issue in the numerical 
noise which is introduced by the (order-by-order) iteration. We explicitly checked that this 
noise was not a big problem (of course FFT-acceleration is quite helpful in reducing the 
number of iterations needed to fix the gauge). In the end the advantages of computing in 
Landau gauge were not overtaken by the care that is due to keep this noise under control. 

We now come to briefly describe how we compute the observablcs of Eq. (1). Trivial 
algebra (i.e. creating external states with quark operators and Wick-contracting to ob- 
tain propagators) leaves us with the task of computing expectation values {i.e. asymptotic 
Langevin time averages) of the quantities 

E M-,],, M-^]^^. (7) 

q;aT 

(where M is the Dirac operator; a and (3 are external polarizations; a and r are other spin 
ind ices; p and q are momentum indices; color degrees of freedom are always suppressed in 
the notation) The index p in the inverse Dirac operator is singled out by placing a 5-like 
source at p in momentum space, with the right polarization and color index (more details in 
the following section). Notice that in this way not only the inverse is to be computed on a 
source (as usual), but one actually squeezes all the information out of the configuration. This 
is the advantage of working directly in momentum space, which is natural in our framework 
(every inversion of M comes as a result of a computation which goes back and forth from 
momentum space; again, see [6] for details). The only measurement which is a bit different 
is that of the conserved vector current 

= 1/2 {^{x) - 1) U,{x) ^{x + //) + ^{x + //) (7;. + 1) ulix) ^{x) ) . (8) 

A little algebra shows that also in this case the measurement can be quite efficient by 
reverting to a convolution product. 

A very important improvement of our statistics comes from exploiting hypercubic symme- 
try: all the measurements connected by a hypercubic symmetry transformation are averaged. 
The fluctuations associated to this average are taken into account for assessing errors. As a 
general rule for the different measurements involved in our calculations, bootstrap was the 
basic tool for the computation of errors. 

To conclude this section about our computational method, we should comment on our 
treatment of the zero modes. Any perturbative expansion of LQCD has to face the problem of 
regularizing the zero modes contribution to the functional integral, since the free propagator 
cannot be inverted in those points. How this applies to NSPT has been discussed in [6], 
where we refer the reader for further details. The most common approach is to remove the 
degrees of freedom associated with the zero modes [14, 1]. Although this prescription is not 
gauge invariant, such contributions are expected to vanish in the inflnite volume limit. We 
should remark that gauge invariant alternatives to this procedure exist. They involve the 
use of twisted boundary conditions [15] or the Schrodinger Functional scheme [5]. While we 
plan to perform computations also in those schemes in the future, in the present work we 
have only considered the prescription in which zero modes are removed. We will come back 
to this issue in section 6. 
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4 Hypercubic symmetric Taylor expansions: the case 
of the quark propagator 



We now proceed to discuss in detail a prototype computation, i.e. the one loop computa- 
tion of the quark field renormalization constant. In practice, we are going to describe how 
we measure the quark propagator. We will thus make clear what we mean by hypercubic 
symmetric Taylor expansions. 

The section is intended as a prototype computation, so let us pin down in general what 
the expected form of the n^^ loop coefficient of a RC is: 

n 

Zn^Cn + J2 log(^)* + F{p) {p = pa) . (9) 

i=l 

We have to look for a finite number (c„), a divergent part which is a function of anomalous 
dimensions 7's and irrelevant pieces, which we can expect compliant to hypercubic symmetry 
and described by a suitable function F. We take the needed anomalous dimensions from the 
literature and we subtract their contribution. In particular, for a one-loop computation we 
simply need to subtract a simple log multiplied by the one-loop anomalous dimension (in 
this section we will completely ignore all the contributions coming from finite-size effects, 
to which we will come back in Section 6). After such a subtraction wc need a convenient 
way to fit the irrelevant pieces given by F{p). The example at hand is both instructive and 
simple: in particular, in Landau gauge the quark field has zero anomalous dimension at one 
loop, so it is simply required that we get rid of F{p) in order to get the constant c„ we are 
interested in. 

We want to compute the two points vertex function (the inverse of the quark propagator) 
for a massless fermion. In the continuum hmit we have: 

r2(/) = s{p')-\ 

On the lattice we define the dimensionless quantity p — pa (in general we use the hat notation 
for dimensionless quantities). Furthermore, we also explicitly write the dependence on the 
coupling (and since we compute in PT we write (3~^ rather than j3) 

ar2{p,mar,/3~^) = aS{p,mar,/3'~^)~^ 

= ip+mw{p) -t.{p,rhcr,P''^) (10) 

where fhw{p) = 0{p^) is the (irrelevant) mass term generated at tree level by the Wilson 
prescription, E(p, rricr, /S"^) is the dimensionless self-energy (which is 0{(3~^)) and rhcr — 
arricr is the critical mass (which is 0{I3'^) as well). Since chiral symmetry is broken by the 
Wilson regularization, also massless fermions generate a mass counterterm. 

The first step is to compute the self energy E from our NSPT simulations. For that we need 
the propagator aSipjihcr, P~^) in momentum space, i.e. 

T 

aS{p, me., = { M-;.^ )=T-'Y1 ^ap;r,p(^) , 

t=l 
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where M^p-^q is the full fermionic matrix. We write explicitly only the spin indices {a and 
f]) and the momentum coordinates (p and q), while color indices are left implicit. The 
symbols () stand for the average over the gauge configurations and the right hand side 
makes exphcit the average over the Monte Carlo history of length T. This is performed as 
described in [6]. Here we only remind that our method - based on a discrctizcd stochastic 
Langevin equation - also involves an extrapolation on the stochastic time discretization. We 
are interested in those elements of the inverse fermionic matrix which appear in the main 
diagonal in momentum space. This is obtained by "sandwiching" the fermionic matrix in a 
(5-like source vector in momentum space: ^^'^\q) — ^aa^pq- The order-by-order inversion is 
then performed as described in [6]. 

Once a5'(p, rhcri is obtained, we average over all the components which are connected 
by hypercubic symmetry transformations. For each given momentum, we numerically (order- 
by-ordcr) invert the 4x4 propagators^. Finally we obtain the self energy E(p, mcr, Z?""^) as 
in Eq. (10). 

Now we turn to the analysis of the self energy that we have obtained as above. It can be 
written as 

E(p, rhcr, = Sc(P, rhcr, P''^) + rricr, P~^) + Sother(P, rhcr: P~^)- (11) 

Sc is the contribution along the (Dirac) identity operator. By this we mean that the trace 
over spin indices: l/4Trspin(5]) = Ec- Similarly, Ey is the contribution along the gamma 
matrices: ^ 

7$l7MTrspin(7M^) = Ey 

Finally, Eothcr includes all contributions along the remaining elements of the Dirac basis. We 
are not interested in such (irrelevant) terms, which are easily projected out. Therefore, we 
will forget about Eothcr in the following. 

Ec contains the contribution to the critical mass, in fact: 

E(0, rhcr, f3~^) = Ec(0, rhcr, P"^) = rhcr = a mar- (12) 

By restoring physical dimensions one can inspect the divergence of the critical mass: 
Ec(0, rhcr, P~^) — "^cr- We will come back to it in the following section. For the moment, 
we concentrate on Ey, which we need to extract the quark field RC. If we make a Taylor 
expansion in powers of a, its most general form up to order O(a^) is: 

±v^iT.^,P, {±P{p,mcr,f3-') +pl±P{p,mcr,p-') + pl±P {p,mcr, P'') + ..•)• (13) 

where the dots stand for higher terms in a. The functions E^ (.) (with i — Q, ...) in turns, 
are the most general combinations of hypercubic-invariant polynomia which contribute to 
the given order. In particular the first term can be written as: 

±Pip, rhcr, r') = + «r + HpI + «f E pIpI + 0{a% (14) 
^Here, we use the fact that the propagator is color diagonal at any order in Perturbation Theory 
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For higher i > 0, there are of course less terms relevant for a given order. In general, 
all the possible covariant polynomial can be found through a character's projection of the 
polynomial representation of the Hypercubic group onto the defining (four dimensional) 
representation of the same group (see for instance [16] for a general reference). 

To gain insight into Eq. (13), remember that in the free case the Sy'' correspond to the 
coefficient of p^*"*"^-* in the Taylor expansion of 2 sin(^). Eq. (13) is what we call a hypercubic- 
invariant Taylor expansion. The term in which we are interested is the leading term af*^ in 
Ey"*, since the other vanish in the continuum limit. In fact, the quark field RC Zq is defined 
as (see Eq. (5); in the following we assume the color average has already been done): 

Trgpin (Xyi/ 

7^p^Sy(p,mcr,/3 i))| 



Zq{iia) = 



In order to explain in details our fit procedure, let us define the auxiliary quantities: 

a{k,p)^^ E ^^P^n(7.Mp,m..,/?-^)) ^^^^ 

{L = Na is the linear size of the the lattice). For a given momentum an average is taken 
over all the M directions /j, such that = ^ k. For instance, consider the momentum 
q — (1, 1,3,2) 2n/N. In this case cr{l,q) has two contributions (M = 2): one from 71 and 
one from 72. Since also the a{k, .) are hypercubic- invariant, they can be averaged accord- 
ingly. Referring to the example above, consider t — (3, 2, 1, 1) 2'!r/N; simmetry requires that 
(t(1, q) = (t(1, t). In practice, they are averaged. Notice that the functions a{k, .) are specific 
linear combinations of the Sy^"*. In practice we fit the data against the functions a{k,.). 
This is nothing but a fit of the constants entering the parametrization of Sy'' in (13). If we 
include 0{a^) terms (as in Eq. (13)), we have to fit 7 unknown constants. To order 0{a^) 
we have 14. We tried different orders up to 0{a^) and check the stability of the result. The 
interesting term is the leading one. Other coefficients have to do with irrelevant effects. 

A nice illustration of the control that we have over our fits is provided by Figure 1. 
There, we plot the functions a{k,p) up to A; = 4, along with the curve Ti^^\p,mcr, I3~^)), 
whose intercept in = is the parameter we are looking for. In Figure 1 we choose to plot 
data versus p^, but one should keep in mind that this is not the only invariant under the 
hypercubic group entering (p, mcr, /S^^)- Our numerical result reproduces very well the 
analytical one. 



5 Results 



In the previous section we saw an example with no anomalous dimension. This is of course 
also the case for finite RC {Zy and Z a) or finite ratios [ZpjZs and Zy IZa). In the following 
wc present our results for these quantities. We computed at every order the relevant expec- 
tations values dictated by Eq. (1). Finally we performed the amputation and the projection 
on the tree level structure. We could thus get the order- by-order expansions of the (9r(pa)'s 
in terms of which RC are defined. One-loop analytical results are well reproduced [17]. 
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Figure 1: The continuum-limit extrapolation of Z^^-* (first loop of the quark field renormal- 
ization constant). We are interested in the intercept at (pa)^ = 0, reached on the lowest 
line, which is the contribution T,y\p,7ficr, f3~^) in Eq. (13). The (blue, violet, azure-blue, 
green) curves represent the functions a{k,p) of Eq. (15) for k = 1,2,3,4. The red curve is 
which is only meant to guide the eye to the intercept at (pa)'^ = 0. The displayed fit 
include up to 0{a^) terms. Stability has been checked with respect to different number of 
terms and intervals. 



In the last subsection our results for the critical mass are presented. 



5.1 The finite ratios Zp/Zs and Zv/Za 



The ratios Zp/Zs and ZyjZA are safely computable at every order. This simply means to 
take (again, order by order) ratios of Or (pa) quantities. The quark field renormalization 
constant present in Eq. (4) drops out in the ratios, together with the divergence that affects 
Zp and Zs separately. In the end, one is left with the same situation we saw in the previous 
section: we simply have to perform at every order hypercubic-invariant Taylor expansions to 
get the continuum- limit coefficients of the expansions. One-loop examples are presented in 
Fig. 2. Fitting a scalar quantity like Zp/Zs is actually easier (there is no direction singled 
out and consequently only one function, to be fitted as a polynomial in the hypercubic 
invariants). 
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Wc could perform many checks on our results. Finite-size effects are well under control, 
as checked by comparing results on 16^^ and 32*^ lattices in the quenched case. In the next 
section we will elaborate on computations for which this is not the case. We also stress that 
we can compute both Z^/Zy and Zv/Za] in the same way, we can compute both Zp/Zs 
and Zs/Zp. Due to the order-by- order nature of the computation, this is not a tautology: 
different ratios come from different (although correlated) combinations of data. We checked 
that to a very good precision the series obtained are inverse of each other. 

Table 1 collects our results for different numbers of flavors. In the case Uf — 2 four-loop 
results are available. As already pointed out, the fact that we were able to go one loop higher 
is due to our better knowledge of the three-loop critical mass in the nj = 2 case. Statistics 
in the cases n/ = 3, 4 is actually poorer. The fact that we could anyway go to three loops is 
a numerical accident: the signals for these ratios are actually very clean. 



Zp/Zs 














o(r') 




o(r') 





- 0.487(1) 


- 1.50(1) 


- 5.72(3) 


n.a. 


2 


- 0.487(1) 


- 1.46(1) 


- 5.35(3) 


- 21.6(3) 


3 


- 0.487(1) 


- 1.43(1) 


- 5.13(3) 


n.a. 


4 


- 0.487(1) 


- 1.40(1) 


- 4.86(3) 


n.a. 


Zv/Za 












O(r^) 


o(r') 




0{(3-') 





- 0.244(1) 


- 0.780(5) 


- 3.02(2) 


n.a. 


2 


- 0.244(1) 


- 0.759(5) 


- 2.83(2) 


- 11.5(2) 


3 


- 0.211(1) 


- 0.711(6) 


- 2.72(2) 


n.a. 


4 


- 0.244(1) 


- 0.732(6) 


- 2.57(2) 


n.a. 



Table 1: The ratios Zp/Zs and Zv/Za for various number of flavor Uf. Four- loop results 
are only available ior Uf — 2. 

Having results for various numbers of flavors one can proceed to flt the n/-dependence. 
Since the polynomial dependence on of every order is flxed, this is another test for our 
results (see Fig. 3). We got 

{Zp/Zsf^^ = -1.50(1) + 0.0249(2) nj {Zp/Zsf''^ = -5.72(3) + 0.151(5) n/ + 0.0159(5) n) 
{Zv/Zj^^^^ = -0.780(5) + 0.0121(1) n/ {Zv/Zj^^^^ = -3.02(2) + 0.073(2) n/ + 0.098(3) n) 

Presented in this (more universal) way the precision of our results appears a bit poorer. 
As expected, results are dominated by quenched contributions. 
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5.2 Zy and Za 



One loop examples of computations of Zy and Za are plotted in Fig. 4. Zy and Za are finite 
quantities by themselves. In our master formula Eq. (4) they are interlaced with log's coming 
from the quark field renormalization constant. The latter can be eliminated in two different 
ways. A first strategy is to cancel Zq directly from the measurements of the propagator. 
Another possibility is to take ratios with the conserved vector current: this is just what we 
did in the case of ZpjZg and Zy/ZA, this time having one of the Z's equal to one. Both 
procedures return consistent results, which are summarized in Table 2, where we present 
results for 71/ = 0, 2 (also in this case, four-loop results are available for n/ = 2). 



Zv 














o(r') 









- 1.044(2) 


- 1.98(3) 


- 6.10(8) 


n.a. 


2 


- 1.011(2) 


- 1.88(3) 


- 5.12(8) 


- 17.0(9) 


Za 












o(r^) 











- 0.800(2) 


- 1.39(3) 


- 4.04(4) 


n.a. 


2 


- 0.800(2) 


- 1.31(3) 


- 3.50(8) 


- 9.8(6) 



Table 2: The finite renormalization constants Zy and for n/ = 0, 2. 

We have just discussed the two different approaches we used to compute Zy and Za- 
In the previous subsection we presented results for the ratio Zy/ZA, which can of course 
as well be computed from the computation of Zy and Za- One can verify that all these 
measurements are very well consistent. Still, they are controlled by different numerical 
noise, so that (for example) a direct computation of the ratio Zy/ZA is viable for all the n/ 
wc took into account, while this is not the case for Zy and Za separately (as already stated, 
statistics for = 3,4 is poorer). In the end, all these procedures difi'cr from each other for 
different ways of fitting irrelevant contributions. Getting rid of irrelevant contributions to 
single out continuum-limit results is a key issue in our approach, and so consistency between 
all these computations is a good test for reliability of our results. 

5.3 A by-product: the critical mass 

Analytical computations of the critical mass are available up to two loops [12]. A three- loop 
computation in the Uf — 2 case was reported by our group in [6]. Here we present three-loops 
result for other Uf and add a four loop result for n/ = 2. Results are collected in Table 3. 
They were obtained from the defining formula of Eq. (12) by fitting irrelevant contributions 
to Sc(P) iTT-cr, P^^)- Also in this case there was no log coming from an anomalous dimension: 
in this case there is a power divergence, in force of which a perturbative result is not to 
be taken as an accurate one. It is nevertheless valuable indeed to maintain massless our 
fermions, i.e. as a counterterm. 
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nf = 


nf = 2 


Hf = 3 


Uf = 4 




- 13.11(6) 


- 11.78(5) 


- 11.02(9) 


- 10.24(9). 


0{(5-') 


n.a 


- 39.6(4) 


n.a 


n.a 



Table 3: Three-loop critical mass for various n/; a four-loop result is available for Uf — 2. 



Also in this case, one can fit a generic Uf result 



= -13.11(6) + 0.62(5) n/ + 0.024(9) 



n 



f 



6 Dealing with anomalous dimensions 

We anticipated that dealing with anomalous dimensions requires some extra care. In order 
to get some insight, we discuss a first example in which an anomalous dimension comes into 
place, i.e. the one-loop computation of Zs- In this case our master formula Eq. (4) reads 




= 1 



in which we explicitly wrote both the constant and the logarithmic contributions to renor- 
malization constants (the only log comes in this case from Zs since one-loop quark-field 
anomalous dimension is zero in Landau gauge). 0^g\p^) is what is actually numerically 
measured. At one-loop order we can solve the previous relation to 



z. 



W-.W=OW(/)-7(^)log(/). (16) 



The message from Eq. (16) is simple: wc will first subtract the logarithmic contribution and 
then proceed to our hypercubic-invariant Taylor expansion. This is plotted in Fig. 5: upper 
data points are O'^^^p'^), lower data points are the subtracted ones. We can see on the left 
of Fig. 5 that by going through this procedure we miss the analytical result. Notice that 
it looks like we were subtracting too much. To be definite, the subtracted data points bend 
quite a lot in the IR region. In the end, this does not come as a surprise: RI'-MOM is 
an infinite-volume scheme, but we are necessarily dealing with finite N (number of lattice 
points) computations. Since for n/ = we have both 32'^ and 16"^ data, we are in a position 
to verify whether this is the real issue. 

Fig. 6 displays our results for 0^\p^) (the equivalent of Eq. (16) for the pscudoscalar 

current), O^^-* {p^) and of the ratio ^li)"^^^ on the two different volumes. While the ratio (in the 

Op (p ) 

middle of the figure) is safe (we have already made this point in the previous Section), quite 
remarkable finite-size effects are manifest for the Of\p^). It is obvious that by performing 
the subtraction of Eq. (16) on the 16"^ data points one misses the analytical result even more 
than in the left of Fig. 5. The picture stays much the same at higher loops. 
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An important caveat is in order. We have already made the point that onr regnlarization 
of zero modes prescribes to remove the degrees of freedom associated to them. This is a 
legitimate procedure in the N ^ oo limit, which in turn means that we can have better and 
better approximations of infinite volume results, but we can not aim at having consistent 
perturbative expansions at finite physical volume. Our aim is to single out the N = oo 
behavior, but this requires to confront finite- A?^ corrections, which we expect to be sizable in 
particular in the IR region. 

One can define L = Na. Let us now write down for the quantity at hand the momentum 
sum I{p, a, L) encoding the lattice Feynman diagram of the conventional Lattice Perturbation 
Theory, with the same ad hoc regnlarization of zero modes (zero momentum removed from 
the sum) . Dimensional analysis suggests the presence of pL — pN effects (this relation holds 
for every value of a). In the spirit of the famous work [18] one can now split a (logarithmically 
divergent) Feynman diagram as in 

I{p, a, L) = /(O, a, L) + (/(p, a, L) - /(O, a, L)) = /(O, a, L) + J(p, a, L). (17) 

We can now manipulate the momentum sums. The divergence is logarithmic so that by 
subtracting /(O, a, L) we make J(p, a, L) UV finite. Therefore it can be computed (with the 
same ad hoc regnlarization of zero modes) in the a — > limit. Although this does not define 
a finite volume perturbative computation, it is a legitimate manipulation of the sum. In 
general, it will be now IR divergent, but this divergence (which is anyway regularized by 
finite L) will be canceled by contributions coming from /(O, a, L), i.e. 

1(0, a, L) = ci + -flog(a/L) + H(a/L) (18) 
J(p, a,L) = C2 + 7 log{pL) + G{pa, a/ L, pL) . 

We point out that /(O, a, L) can not contain pL effects: these should be looked for in 
J{p,a,L). Therefore one can look for pL — pN effects in G{pa,a/ L,pL) — > G{pL). In 
order to obtain this quantity, we just computed the relevant graph in the formal continuum 
limit of our sum J(p, a, L) (a — > with L = Na fixed), with the same ad hoc regnlarization 
of zero modes. We call this contribution a tamed-log, since it is supposed to resemble the 
expected log, but with pL = pN effects on top of it. We find that this function indeed 
approaches a log for p >> 1. Fig. 5 displays our results once one subtracts this tamed-log. 
As a matter of fact, if one stays away from deep IR the subtracted data points on the left 
and on the right of Fig. 5 are much the same. We stress that we are not saying that the 
finite effects we have just elaborated on are the only ones. By inspection, they appear to 
be the relevant ones, as it confirmed by the fact that A^ = 32 and A^ = 16 now return the 
same results. 

The situation is more complicate at higher loops. We will devote to it a separate paper 
[9], in which we will explicitly gain informations from different lattice sizes. 



7 Resumming the series 

We now go back to the expansions of subsection 5.1 and 5.2 and try to resum them to 
obtain the finite RC's. Giving results and errors on top of them requires the estimation of 
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truncation errors. Wc will in the following adopt the strategy of BPT. We stress from the 
very beginning that our real goal is the estimation of convergence properties of the series. It 
is only in force of the sufficiently high order of the expansions that one can hope to really gain 
insight. One should nevertheless be ready to accept that every statement on convergence 
will be decided on a strict case-by-case pohcy. 

The different coupling constants we will use are all obtained in terms of the basic plaquette 
P. Let us define 

xq = p-^ xi = ^2 = "2 log(^) ^3 = -p- ■ (19) 

X2 and ,T3 arc quite popular as boosted couplings. The reason why wc also define Xi will be 
clear in a moment. Obtaining the expansions in Xi once the expansions in xq are known is a 
textbook exercise, given the definitions in Eq. (19). One needs the expansion of the plaquette: 
analytical results [19] are only known to a given order, but our simulations always provide 
also the expansion of P. 

We resum the series at /3 = 5.8, n/ = 2. This makes possible a comparison with the 
non-perturbative results of [8]. 

Fig. 7 displays the resummation of Zp/Zs and Zs/Zp in the four different couplings. One 
can inspect from the very beginning the impact of a basic property of BPT which is often 
underestimated: all the couplings are equal at tree level, which means that all the expansion 
are equal at leading order. One loop BPT amounts to sitting on a straight line, whose slope 
is dictated by the one-loop coefficient. Only at higher loops we can gain some insight on 
convergence properties. There is actually a variety of convergence patterns (taking also Xi 
into account is helpful with this respect). 

In particular, one can check the following: 

• Within a fixed definition of the coupling, convergence is of course better and better as 
the order increases. As common wisdom suggests, convergence in the bare coupling is 
not so brilliant and in general quite different convergence patterns are manifest; they 
appear quite satisfactory for X2 and 2:3. In particular, for the case of X2 expansion 
we plot in Fig. 8 the deviations A*^"'^ defined as the differences between resummation 
at order n and resummation at order n — 1 . The good scaling should not be taken 
too seriously (this is largely a numerical accident) . Still, this is signaling a reasonable 
convergence pattern. 

• As the order increases, expansions in different couplings get closer to each other, as 
expected; in particular expansions in X2 and X3 are quite close to each other. 

• The rcsummed results for Zp/Zs and Zs/Zp in the X2 and ,T3 couplings are the inverse 
of each other to a reasonable accuracy. This is also a good indication. 

Convergence properties of the expansions in the X2 and 2:3 couplings are good enough 
to extract a result. We notice that if one adds to the result at a given order the deviation 
from the immediately lower order, one always ends up at the same result (as a matter of 
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fact a popular way to pin down a truncation error is just taken from deviations which we 
previously called A„). We thus quote ZpjZg = 0.77(1). 

We have already made the point that to assess convergence properties one should adopt 
a case by case strategy. This can be clearly seen when we proceed to resum Zj^ and Zy. 
Before doing that, we give a trivial example of what a blind application of the idea of BPT 
can result in. In Fig. 9 we exaggerate the boosting of the coupling, by taking into account 
other coupling Xa = (a > 1). As one can see, convergence properties are completely 
jeopardized. 

In Fig. 10 we plot the resummation of Zy and Za (again, at P = 5.8, n/ = 2). As one can 
see, this time convergence properties of the expansion in the bare coupling are not so bad. 
Consequently, one is already at risk of overshooting at one loop BPT and the expansions in 
X2 and ^3 oscillate. Our final estimates are Za — 0.79(1) and Zy = 0.70(1). 

Our resummed results are quite consistent with [8]. A bigger deviation is seen on the 
values of Za and Zy. To our understanding this could be mainly imputed to the indetermi- 
nation coming from the chiral extrapolation. 



8 Conclusions 

We presented high-order computation of renormalization constants for Lattice QCD. Finite- 
size effects are well under control for the quantities we considered. There is no extrapolation 
involved in staying at the chiral limit in which renormalization conditions are imposed. The 
continuum limit extraction is achieved in a clean way. Truncation errors can be well assessed 
by a judicious use of BPT. Thus, the main message of this paper is that high precision 
perturbative computations of lattice QCD renormalization constants are feasible and should 
not be regarded necessarily as a second choice. 

Further work will follow, both to complete the job for logarithmically divergent quantities 
and to take into account different actions (in particular different fermionic regularizations) . 
This is not expected to imply any change in strategy and the implementation is mainly a 
matter of programming. In particular, work has already started to extend results to Clover 
fermions and to other gauge actions. 
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Figure 2: Computation to one loop of finite ratios of renormalization constants: Zp/Zs (left) 
and Za/Zy (right). Data points taken into account in these particular fits are enclosed in 
circles (left) or joined by solid lines (right; see caption of Fig. 1). 
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Figure 4: Computation to one loop of finite renormalization constants: Zy (left) and Za 
(right). Same notations as in Fig. 2. 
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Figure 5: Computation of one loop renormalization constant for the scalar current. With 
respect to Eq. (16), upper points are the unsubtracted Of^ij?), while lower (circled crosses) 
stand for the subtracted O^^^(p^) — 7^^-* log(]5^). Analytic result is marked with a darker 
symbol. On the left: no correction for finite volume. On the right: finite-volume tamed-log 
taken into account. 
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Figure 6: Computations of 0^\p'^) (the equivalent of Eq. (16) for the pseudoscalar current) 
(top) and 0^^\p'^) (bottom) on 32^ (circles) and 16^ (diamonds). In the middle the ratio 



-, which appears safe with respect to finite-size effects. 
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Figure 7: Resummations of Zp/Zs (left) and Zs/Zp (right) for n/ = 2 at /5 = 5.8 to one 
(circles), two (squares), three (diamonds) and four (triangles) loops (the last is the only one 
which has a sizable error). We show resummations for different couplings: on the x-axis, the 
(different) values of the different couplings. From the left: xq, xi, X2, X3 (xq is see text 
for the definitions of the other couplings) . 
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Figure 8: The scaling of deviations of different order truncations for the quantity Zp/Zs for 
the X2 coupling. 
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Figure 9: Resummations of Zp/Zs (left) for Uf = 2 dX j3 = 5.8. The same as in Fig. 7, but 
this time exaggerating the boosting of the couphngs. 
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Figure 10: Resummations of Z^ (left) and Zy (right) for n/ = 2 at /3 = 5.8 to one (circles), 
two (squares), three (diamonds) and four (triangles) loops (the last is the only one which has 
a sizable error). We show resummations for different couplings: on the x-axis, the (different) 
values of the different couplings. From the left: xq, (see text for the definitions of the 

couplings) . 
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